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Three coarse-grained molecular dynamics (MD) 
models are investigated with the aim of developing 
and analysing multi-scale methods which use MD 
simulations in parts of the computational domain and 
(less detailed) Brownian dynamics (BD) simulations 
in the remainder of the domain. The first MD model 
is formulated in one spatial dimension. It is based 
on elastic collisions of heavy molecules (e.g. proteins) 
with light point particles (e.g. water molecules). Two 
three-dimensional MD models are then investigated. 
The obtained results are applied to a simplified 
model of protein binding to receptors on the cellular 
membrane. It is shown that modern BD simulators 
of intracellular processes can be used in the bulk 
and accurately coupled with a (more detailed) MD 
model of protein binding which is used close to 
the membrane. 



1. Introduction 

Brownian dynamics (BD) simulations have been used for 
the modelling of a number of spatio-temporal processes 
in cellular and molecular biology in recent years, 
including models of intracellular calcium dynamics 
[1], the MAPK pathway [2] and signal trasduction in 
Escherichia coli chemotaxis [3]. In these applications, 
trajectories and interactions between key biomolecules 
(e.g. proteins) are calculated using BD methods, while 
other components of the system (e.g. solvent molecules), 
which are of no special interest to a modeller, are not 
explicitly included in the simulation, but contribute to the 
dynamics of Brownian particles collectively as a random 
force. This reduces the dimensionality of the problem, 
making BD less computationally intensive than the 
corresponding molecular dynamics (MD) simulations. 
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Denoting the position of a Brownian particle by X = [Xi,X2,X3] and its diffusion constant by 
D, a simple model of Brownian motion is given by the (overdamped) Langevin equation 

dX z = VlDdWf, i = 1,2,3, (1.1) 

where Wj, i — 1,2,3, are three independent Wiener processes [4]. BD approaches which are 
based on (1.1) have been implemented in a number of software packages designed for spatio- 
temporal modelling in systems biology, including Smoldyn [5], MCell [6], Green's function 
reaction dynamics [7] and the first-passage kinetic Monte Carlo method [8]. The software package 
Smoldyn discretizes (1.1) using a fixed time step At, i.e. it computes the time evolution of the 
position X = X(t) of each molecule by 

Xi(t + At) = Xi(t) + VlDAt % it i = 1, 2, 3, (1.2) 

where [^1,^2/^3] is a vector of normally distributed random numbers with zero mean and unit 
variance. A different BD approach is implemented in the Green's function reaction dynamics [2] 
which evolves time using a variable time step. It approximately computes the time when the next 
reactive event happens. This means that trajectories of molecules which are not surrounded by 
other reactants can be simulated over longer time steps. 

Although the BD models are becoming a popular choice for stochastic modelling of 
intracellular spatio-temporal processes, several difficulties prevent the use of BD for some 
systems. First of all, detailed BD models are often more computationally intensive than coarser 
spatio-temporal models which are written for concentrations of biochemical species. In some 
applications (e.g. intracellular calcium dynamics [1] or actin dynamics in filopodia [9]), individual 
trajectories (computed by BD) are important only in certain parts of the computational domain, 
while in the remainder of the domain a coarser, less detailed, method can be used. In these 
applications, the computational intensity of BD simulations can be decreased by using multi- 
scale methods which efficiently and accurately combine models with a different level of detail in 
different parts of the computational domain [10,11]. 

Another difficulty of BD simulations in cell and molecular biology is that detailed BD models 
require more parameters than coarser (macroscopic) models. In some studies, macroscopic 
parameters are used to infer BD parameters [5,12]. For example, knowing the macroscopic 
reaction rate k of a bimolecular reaction A + B -> C and diffusion constants of reactants, one 
can calculate a (microscopic) reaction radius of BD simulations which gives the corresponding 
macroscopic parameters in the limit of many particles. In the classical Smoluchowski limit 
[13], a bimolecular reaction occurs whenever the distance of reactants is less than the 
reaction radius 

* = 4*(Da + Db)' (L3) 

where Da (resp. Db) is the diffusion constant of reactant A (resp. B). Although this approach is 
commonly applied in stochastic reaction-diffusion models, it is not the most satisfactory, because 
different microscopic models can lead to the same macroscopic process and parameters [12,14]. 
For example, the simplest Smoluchowski model (1.3) assumes that all collisions are reactive but, 
in reality, many non-reactive collisions of molecules happen before a reactive collision occurs. 
Therefore, some algorithms postulate that molecules only react with a certain rate (probability) 
when the distance between reactants is less than a modified reaction radius (which is larger 
than q). Other methods discretize the Langevin equation with time step At and substitute the 
Smoluchowski formula (1.3) (which is valid for an infinitely small time step) by a tabulated 
function computed numerically [5]. However, all of these approaches are verified by considering 
the macroscopic limit (of many reactants) and showing that the reaction occurs with the given 
rate k in this limit. 

A different approach to parametrize BD models is to use a more detailed description written 
in terms of MD. In this paper, we investigate connections between BD and MD models with the 



aim of developing and analysing multi-scale methods which couple BD and MD simulations. 
We consider a (computationally intensive) MD simulation in a domain Q which is either 
one or three dimensional, i.e. Q C K or Q c M 3 . Our main goal is to design and analyse multi- 
scale methods which can compute spatio-temporal statistics with an MD level of detail in the 
subdomain C & . We define 



where and Q are open sets, the (open) set Qq is the complement of and I is the shared 
interface (boundary) between and Qq. In the multi-scale set-up (1.4), we use a detailed MD 
model in fi^ and a coarser BD model in Qq. 

In this paper, we focus on a simple MD approach which is introduced in §2 and §3. A few 
(heavy) particles with mass M and radius R are coupled with a large number of light point 
particles with masses m<^M. The collisions of particles are without friction, which means that 
post-collision velocities can be computed using the conservation of momentum and energy 
[15,16]. We will introduce and study three MD models which make use of elastic collisions. They 
will be denoted as MD models [A], [B] and [C] in what follows. More complicated MD approaches 
are discussed in §7. 

The first MD model [A] is introduced in §2. It is a one-dimensional MD model where all 
particles move along the real line. In particular, the radius R does not have to be considered, 
because it has no influence on the dynamics of large particles. In one dimension, heat bath 
particles cannot pass each other, which makes the MD model [A] different from three-dimensional 
models in §3 where heat bath particles (points) do not interact with each other. 

In §3, we introduce two three-dimensional models, denoted [B] and [C], where the non-zero 
radius R is one of the key parameters. To make one- and three-dimensional models comparable, 
we keep R fixed in the three-dimensional model and we study the behaviour of all MD models in 
the limit M/m -> oo. This limit can be achieved in many different ways. For example, we can keep 
m fixed and pass M -> oo, or we can keep M fixed and pass m -> 0. In what follows, we define the 
parameter 



This parameter is dimensionless, even if we assume that M and m have physical units of mass. 
However, in this paper, all parameters are considered dimensionless for simplicity. We are 
interested in the limit \x -> oo. 

All three models [A], [B] and [C] converge in appropriate limits to the Brownian motion of 
large particles given by (1.1). One can also show that these models converge to the Langevin 
description [15-17] 



where [Xi,X2,X3] is the position of a diffusing molecule, [V\, V^,!^] is its velocity, D is the 
diffusion coefficient and y is the friction coefficient. This description can be further reduced 
to (1.1) in the overdamped limit y -> oo. We overview the results which relate MD models [A], 
[B] and [C] with Brownian motion in §2 and §3. 

Both (1.1) or (1.6)-(1.7) reduce the dimensionality of the problem, making BD less 
computationally intensive than the corresponding MD simulations. In §4 and §5, we study how 
MD models [A], [B] and [C] can be used in one part J2d of the computational domain Q and the 
BD models (1.1) or (1.6)-(1.7) in the remainder Qq, making use of the notation (1.4). We apply our 
findings to a simplified model of protein binding to receptors in §6. We conclude with discussing 
our results in §7. 



Q c = Q \ J2 D and I = d£2 D n dQ Cf 



(1.4) 



dX; = Vi dt, 



(1.6) 




(1.7) 



2. One-dimensional molecular dynamics model [A] 

The MD model [A] is described in terms of positions x l and velocities v l , i = 1,2,3, . . ., of heat 
bath particles, and positions X 1 and velocities V 1 , i = l,2, . . . ,N, of heavy particles of mass 
M^>m, where m is the mass of a heat bath particle [15]. In our computer implementations, we 
will consider a finite number of heat bath particles. However, we formulate the MD model in 
terms of (countably) infinitely many heat bath particles which are initially distributed along the 
real line according to the Poisson distribution with density 



where \x is given by (1.5), and D and y are positive constants. This means that the probability that 
there are j particles in a subinterval [a, b] cR, a <b, is equal to 



- exp [ - X^b - a)], 



where (b — a) is the length of the interval [a, b\. Initial velocities of heat bath particles are given by 
the normal distribution 

fM = — U= ex P ( " ^2 V where °> = y/fa + WY' (2-2) 

Let us consider a model with a single heavy particle, i.e. N = 1. Then its location X 1 and velocity 
V 1 will be denoted as X and V to simplify our notation. Whenever the heavy particle collides with 
the light particle with velocity v l , their velocities are updated using the conservation of mass and 
momentum by 

V=^V+T^" i (2-3) 
M + m M + m 

and , K r. , , 

; m-M : 2M T7 

v = Tl v + Tl V > 2 - 4 

M+m M+m 

where tildes denote post-collision velocities. Using (1.5), the equations (2.3) and (2.4) can be 

rewritten as . 1 

V=^iV+^-v i and = t—^v* + -=!±-V. (2.5) 
li + l fi + 1 /x + 1 /x + 1 

The following result can be shown for the above MD model [A]. 

Lemma 2.1. Let y > 0 and D > 0. Let us consider the heavy particle of mass M with initial position 
X /x (0) = Xq and initial velocity V /ji (0) = Vo which is subject to elastic collisions (2.5) with heat bath 
particles of mass m whose initial positions and velocities are distributed according to (2.1) and (2.2). Then 
X M and V M converges (as \i -> oo) in distribution to the solution X and V of equations 

dX=Vdt and dV = -yV + y V2DdW, (2.6) 

where X(0) = Xq and V(0) = Vq. That is, X and V solve the one-dimensional version of equations (1.6) 
and (1.7). 

Proof. This lemma can be proved using the main theorem in Holley [15], where it is shown that 
a similar process converges to the Ornstein-Uhlenbeck process (1.7) for velocities. Although our 
function/^ does not satisfy all assumptions of the main theorem of Holley [15], a simple rescaling 
of our parameters leads to a process which is covered by Holley's theorem. In §4 of this paper, we 
also rederive this result as one of the consequences of multi-scale analysis (see (4.11)). ■ 

As the goal of this paper is to study the behaviour of computational algorithms, we formulate 
the MD model [A] in a finite domain [— L, L], i.e. we consider a finite number n = n(t) of heat bath 
particles which are at positions x l e [-L, L] with velocities v l e (—oo, oo), i = 1, 2, . . . , n. We want to 
formulate boundary conditions of our problem so that the spatio-temporal statistics in [— L, L] are 
equivalent to spatio-temporal statistics of the original unbounded process. The following lemma 
will be useful for designing appropriate boundary conditions. 



Lemma 2.2. Let b e R and At>0. Let us assume that heat bath particles are distributed according 
to the Poisson distribution with density (2.1) in the interval (—oo,b). Their initial velocities are given 
according to (2.2) and there are no particles in the interval (b, oo) at time t = 0. Then the average number 
of particles in the interval (b, oo) at time t = At is 

y(M + l)A£ 

8 { ' } 

The positions x and velocities v of these particles are distributed according to 

H(b-x+vAt)X^(v), (2.8) 

where H(-) is the Heaviside step function. In particular, the positions of the particles at point x e (b, oo) are 
distributed at time t = At according to 

X I % — b \ 

q(x; At, b) = -^- erfc , for x e (b, oo), (2.9) 

2 y<7 M AfV2/ 

where erfc(z) = exp(— s 2 ) ds is the complementary error function. 

Proof. Particles which are at point x e (b, oo) at time t = At were previously at point x — v At at 
time t = 0. In particular, there will be non-zero heat bath particles with velocity v at point x at 
time t = At provided that x — v At < b, which implies (2.8). Consequently, the density of particles 
which are at point x e (b, oo) at time t = At is 



Q(x;At, b) = 



H(b - x + vAt)X^(v) dv = 

X t 



X^f^dv 

(x-b)/At 



7== exp I - \dv = ^- erfc ( — — ^— ] . 



Thus, we proved (2.9). Integrating this formula over x in interval (b, oo), we obtain the average 
number of particles which are in the interval (b, oo) at time t = At 



> x 
Q(x;At,b) dx=-^ 



coo 



erfcl^^—Jdz^^: 
0 Xa^AtVi/ V2tt 



XnCTnAt 



Substituting (2.1) for A M and (2.2) for o^, we obtain (2.7). ■ 

We use lemmas 2.1 and 2.2 to design a computational test for multi-scale methods. As the 
number n = n(t) of heat bath particles in [— L,L] is much larger than the number N of large 
particles, we will focus on models of a single large particle, i.e. N = 1, which is described by 
its position X and velocity V. We choose a small time step At. One iteration of the MD algorithm 
is presented in table 1. We first compute the positions of all particles at time t + At in step [Al] by 
assuming that particles do not interact. Then we use (2.5) to incorporate collisions in step [A2]. 
As all heat bath particles have the same mass, the collisions between them result in exchange 
of colliding particles' positions and velocities. In particular, step [A2] can be implemented by 
sorting the heat bath particles during every iteration. If the large particle collides with a heat bath 
particle, then their positions at time t + At will be obtained in step [A2] by calculating the exact 
time t c of their collision, i.e. t c e[t,t + At]. Then we use their pre-collision velocities in the time 
interval [t, t c ] and post-collision velocities in the time interval [t c , t + At] to compute X(t + At) 
and x l (t + At). All heat bath particles which left the domain [— L, L] are removed in step [A3]. 

New heat bath particles are introduced in steps [A4] and [A5]. We assume that At is chosen so 
small that (2.7) is much smaller than 1. Then (2.7) can be interpreted as a probability of introducing 
one particle from the left (resp. right) during one time step. Using lemma 2.2, the new particle 
will be introduced at the (left boundary) position which is sampled according to the probability 
distribution proportional to q(x; At, —L) in step [A4]. To sample from this probability distribution, 
we scale and shift a random number sampled from the complementary error function distribution 



Table 1. One iteration of the computer implementation of MD model [A]. 



6 



b 
(3 



[A1] Compute 'free-flight positions' of heat bath particles and the large particle at time t + At by 

x\t + AO = x\t) + v\t) At and X(t + At) = X(t) + V(t) At. 

[A2] Compute post-collision velocities by (2.5) for every pair of particles which collided. Compute their post- 

collision positions x\t + At) and X(t + At) by updating their 'free-flight positions' x l (t + At) and 
X(t + At). 

[A3] Terminate trajectories of heat bath particles which left the domain [-L, L]. Update n accordingly. 

[A4] Generate a random number r\ uniformly distributed in (0, 1). 

If ri < y(fi + l)A£/8, then increase nbyl, and introduce a new heat bath particle at a position sampled 
according to the probability distribution proportional to q(x; At, —L). Its velocity is sampled according to 
the probability distribution proportional to H(—L — x + vAt)f^(v). 

[A5] Generate a random number r 2 uniformly distributed in (0, 1). 

Ifr 2 < + l)A£/8,then increase n by 1, and introduce a new heat bath particle at position x n (t + 
At) with velocity v n (t + At) which are sampled according to probability distributions (2.11) and (2.12). 

[A6] Continue with step [A1] using time t = t + At. : g 

; uu 

Table 2. The acceptance-rejection algorithm which is used to sample random numbers which are distributed according to the 
probability distribution n erfc(z), where z e (0, oo). In our simulations, we use constants ^ and a 2 given by (2.10). 

— Generate a random number £i uniformly distributed in (0,1). 

— Compute exponentially distributed random number £ 2 by £2 = — 01 log(£i). 

— Generate a random number £3 uniformly distributed in (0,1). 

— If £1 £3 < #2 erfcfe), ^en choose £ 2 as a sample from the probability distribution n erfc(z). Otherwise, repeat the 
algorithm. 

7i erfc(z), where z e (0, 00). An acceptance-rejection algorithm for sampling random numbers from 
71 erfc(z) is given in table 2. We use it with the constants a\ and given by 

01= 0.532 and a 2 = 0.814. (2.10) 

The values of constants a\ and a 2 were numerically estimated to maximize the total acceptance 
probability 02/(^1 V^ 7 ) °f the acceptance-rejection algorithm in table 2. Using (2.10), we accept 86% 
of proposed numbers £2- 

A particle introduced close to the right boundary in step [A5] will have its position sampled 
according to the probability distribution 

Ci erfc ( L ~* |, for x g (-00, L), (2.11) 

where C\ is a normalization constant. The probability distribution (2.11) is proportional to 
q(—x; At, —L) and can be justified using the same argument as lemma 2.2. To sample from the 
probability distribution (2.11), we again use the acceptance-rejection algorithm in table 2 with 
parameters a\ and a 2 given by (2.10). In step [A5], we also sample the velocity v e R of the new 
particle using the truncated Gaussian distribution 

C 2 H(x-vAt-L)f fl (v), (2.12) 

where C 2 is a normalization constant. To sample random numbers according to the truncated 
normal distributions in steps [A4] and [A5], we use an acceptance-rejection algorithm which is 
derived as proposition 2.3 in [18]. 
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Figure 1. (a) Thirty illustrative trajectories of the heavy particle computed by the MD algorithm [A1]-[A6]. (b) The mean 
squared displacement computed by 10 3 realizations of the algorithm [A1]-[A6] (solid line). The MD results are compared with 
BD results: equation (2.13) (dashed line), Vwi (dot-dashed line) and equation (2.14) (dotted line). We use /x = 10 3 , y = 10, 
D = 1, At = 1(T 7 , L = 20, X(0) = 0 and l/(0) = 0. (Online version in colour.) 



In figure 1, we present illustrative results computed by the algorithm [A1]-[A6]. We use 
\x = 10 3 , y = 10 and D = 1. We initialize the position and velocity of the heavy particle as X(0) = 0 
and V(0) = 0 and we use the algorithm [A1]-[A6] with time step At = 10 -7 in the interval 
[— L,L] where L = 20. In figure la, we present 30 illustrative trajectories of the heavy particle X(t) 
computed for t e [0, 10]. The mean squared displacement given by the MD model [A] is plotted in 
figure lb as the solid line. To illustrate the limiting result in lemma 2.1, we also plot the mean 
squared displacement corresponding to the limiting solution X of (2.6). It can be analytically 
computed as 

MX® - Xm^llDt - ™ + 4Dex Pt-^] _ Dexp[-2^ 

v v y Y Y 

where E[«] denotes the expected value. It is plotted as the dashed line in figure lb. We also plot 
the mean squared displacement corresponding to the overdamped limit (1.1), i.e. \l7Dt, as the 
dot-dashed line in figure lb. If we neglect the exponential terms in (2.13), we obtain 



7e[(X(0-x(0))2] 




(2.14) 



This approximation is plotted in figure lb as the dotted line. We will use (2.14) later in §6 to couple 
the overdamped BD model (1.1) with MD simulations. 



3. Three-dimensional molecular dynamics models [B] and [C] 

MD models [B] and [C] are three-dimensional generalizations of the MD model [A]. They 
are described in terms of positions x l and velocities v l , i = 1, 2, 3, . . ., of heat bath particles, 
and positions Xj, = [X^ ;1 , X^ ;2 , X^ ;3 ] and velocities Vj, = [V* v V^ ;2 , V^ ;3 ], i = 1, 2, . . . , N, of heavy 
particles of mass My>>m, where m is the mass of a heat bath particle [16]. We again define \x 
by (1.5). We will denote by R the radius of a heavy particle. 

MD models [B] and [C] are both based on elastic collisions of heavy molecules (balls with mass 
M and radius R) with point bath particles with masses m. As the collisions are without friction, 
conservation of momentum and energy then yields the following generalization of formulae (2.5) 



for post-collision velocities [16] 

^ = r<]ii + ^[v^ + -^[v;]\ (3.i) 

v7 = [v/]ll + + ^ T [V| i ] ± , (3.2) 

/X + 1 [l + 1 

where is the velocity of the heat bath molecule which collided with the z'th heavy molecule, 
tildes denote post-collision velocities, superscripts _L denote projections of velocities on the line 
through the centre of the molecule and the collision point on its surface, and superscripts || denote 
tangential components. 



(a) Molecular dynamics model [B] 

MD model [B] will use the normal distribution for velocities of heat bath particles. The following 
lemma generalizes lemma 2.1 to the three-dimensional MD model [B]. 

Lemma 3.1. Let y > 0, D > 0 and R>0. Let us consider the MD model [B] where heat bath particles 
are distributed according to the Poisson distribution with density 



* 8R 2 V IttD ' 
Let the velocities of heat bath particles be distributed according to 



3 ^ + (3.3) 



^ (V)= ^)372 eX P 



9 o 9 ' 

v^ + v^ + V3 

H 



where = y/(ji + l)Dy (3.4) 



and v = V2, ^3]. We will consider one heavy molecule in such a heat bath, i.e. N = 1. Then the position 
and velocity of the heavy molecule, X M and V M , converge (in the sense of distribution) to the solution of 
(1.6) and (1.7) in the limit \i -> 00. 

Proof. The MD model [B] and heat bath distributions (3.3) and (3.4) satisfy the assumptions of 
theorem 2.1 in Diirr et al. [16]. Their theorem expresses the limiting equation of a process with 
given A M and/ M (v) in terms of moments of/ M . These moments can be analytically evaluated to 
verify the statement of lemma 3.1. We will also rederive this result in §5 as a consequence of the 
analysis of multi-scale methods. ■ 

Lemma 3.1 can be viewed as a different formulation of theorem 2.1 in [16]. They were interested 
in the limit m -> 0, which is equivalent to \x -> 00. Considering the scaling m^^vm 1 / 2 ) of the 
velocity distribution of heat bath particles (with density scaled as X/m 1 / 1 ), they derived formulae 
for y and D in terms of moments off and k. To formulate lemma 3.1, we inverted their results by 
deriving the appropriate distributions (3.3) and (3.4) which lead to the limiting BD model with a 
given D and y . 



(b) Molecular dynamics model [C] 

In lemma 3.1, we used the normal distribution for velocities (3.4). Another option is to use heat 
bath particles with fixed speed as is done in the following lemma 3.2. We denote the resulting MD 
model as the MD model [C]. 

Lemma 3.2. Let y > 0, D > 0 and R>0. Let us consider the MD model [C] where heat bath particles 
are distributed according to the Poisson distribution with density 



3 I (H + 1)Y 
</x 8ttR 2 V D 



(3.5) 



Let the velocities of heat bath particles be distributed according to 

= ^2 8 \s/ v l + v 2 + v 3 ~ ' ^herecr fJi =2y/(fi + l)Dy (3.6) 

and 8 is the Dirac delta function. Let us consider one heavy molecule in this heat bath at position X M with 
velocity V M . Then X M and V M converge (in the sense of distribution) to the solution of (1.6) and (1.7) in 
the limit \i -> oo. 

Lemma 3.2 can again be proved using theorem 2.1 in [16] that is applicable to any spherically 
symmetric velocity distribution which has at least four finite moments. 

(c) Boundary conditions for molecular dynamics models [B] and [C] 

Next, we generalize lemma 2.2 to the three-dimensional case. This will help us to specify 
boundary conditions for simulations which use the MD models [B] and [C] in finite domains. 

Lemma 3.3. Let b eR and At > 0. Let us assume that heat bath particles are distributed according to the 
Poisson distribution with density A M in the half space (— oo, b) x R 2 ; their initial velocities are distributed 
according tof^v) and there are no particles in the half space (b, oo) x R 2 at time t = 0. Let us assume that 
A M andf^iy) are either given by (3.3) and (3.4) (MD model [B]) or given by (3.5) and (3.6) (MD model 
[C]). 

Then the positions x and velocities v of heat bath particles in the half space (b, oo) x R 2 are distributed 
at time t = At according to 

Hty-xx + viAt)^^), (3.7) 
and the average number of particles in the semi-infinite cuboid (b, oo) x (0, l) 2 at time t = At is 

3y(n + l)At 

16nR2 ■ (3 - 8) 
Proof. Formula (3.7) is a generalization of formula (2.8) in lemma 2.2 and can be justified using 
the same arguments. To prove (3.8), we will distinguish two cases. 

First, let us consider that and / M (v) are given by (3.3) and (3.4). Integrating (3.7) over 
positions and velocities (see the proof of lemma 2.2), we conclude that the average number of 
particles in the semi-infinite cuboid (b, oo) x (0, l) 2 at time t = At is equal to 

3y(/x + l)Af 



^■fxer^At 



16ttR 2 



which is formula (3.8). 

Next, let us consider that A M and/ M (v) are given by (3.5) and (3.6). Integrating (3.7) with respect 
to v, we get the density of particles at x e (b, oo) x R 2 at time t — At, 



Q(x;At,l 



H(b-x 1 + v 1 At)k IJ f ll (v) dv 
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fx 



X " (a 



(xi-b)/At 



v\ + v\ + v\ - cr^ dv2 dv^ ^ dvi 



At 



(3.9) 



where (•)+ denotes a positive part. Integrating this formula over x in the semi-infinite cuboid 
(b, oo) x (0, l) 2 , we obtain 



q(x; At,b) dx = 



(fc,oo)x(0,l) 2 



Substituting (3.5) for A M and (3.6) for o^, we obtain (3.8). 
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Figure 2. (a-c) Schematic of one-dimensional multi-scale set-up (1.4). {d) Schematic of multi-scale set-up (1.4) (in two 
dimensions). (Online version in colour.) 

Lemma 3.3 can be used to specify boundary conditions for simulations of the MD models [B] 
and [C] in finite domains as we did for the one-dimensional case in lemma 2.2. In §5, we will use 
lemma 3.3 to develop and analyse multi-scale approaches which can efficiently and accurately 
compute results with an MD level of detail in a (relatively small) sub domain c £2 by using 
coarser BD simulations in the remainder. The geometry of the desired multi-scale method is 
formulated using (1.4) where an MD model is used in a coarser BD model is used in £2q and 
these models are coupled across the interface L The situation is schematically shown in figure 2d, 
which presents a two-dimensional version of our multi-scale set-up. Here, point particles describe 
heat bath molecules which are used in Q^. Large biomolecules of interest are denoted as grey 
circles. They are simulated using BD in Qq. The solid line denotes interface I. 

The schematic in figure 2d is presented in two spatial dimensions to better visualize the 
problem geometry. MD models [B] and [C] are formulated in a three-dimensional physical space. 
In the three-dimensional version of figure 2d, the cloud of point particles would cover the grey 
ball. To get some insights into this multi-scale problem, we start with the one-dimensional MD 
model [A]. 

4. From one-dimensional molecular dynamics model [A] to Brownian dynamics 

In the case of one-dimensional MD model [A], the situation is schematically shown in figure 2a-c, 
where we consider only one large (heavy) particle, i.e. N = 1. The large particle can either be in 
Qq (figure 2a), or be in fi^ (figure 2c) or crossing the boundary as is shown in figure 2b. Our 
geometry is given by (1.4), where 

Q = (-L,L), £2 D = (-L,0), J2 C = (0,L) and I = {0}. 

The large particle covers the interval (X(t) — R,X(t) + R). Let us consider that the large particle 
intersects the interface I, as is shown in figure 2b. Then I c (X(t) — R, X(t) + R), which is equivalent 
to X(t) e (—R,R). The heat bath particles are simulated in using the MD model [A]. Let us 
choose At so small that the probability of two collisions happening in the time interval (t, t + At) 
is negligible. As we do not explicitly simulate heat bath particles in Qq, we will consider an 
additional correction of the velocity of the heavy particle in the form 

V(t + At) = V(t + At) + a(V(t))At + 0(V(f))VAf f , (4.1) 

where V(t + At) is the post-collision velocity of the heavy particle at time t + At, which only 
takes into account collisions with the heat bath particles from the left. It is either equal to V(t) or 
computed by (2.5) if a collision with a heat bath particle occurred in Q^. Equation (4.1) is adding 
both the drift term a(V(t))At and the noise term p{V{t))\[~Ai% , where £ is a normally distributed 
random number with zero mean and unit variance. The drift and noise terms implicitly take 



into account collisions at the right boundary (X(t) + R) of the heavy particle. Passing At -> 0, we 
observe that the contributions of the collisions at the right boundary are given by the Ito stochastic 
differential equation 

dV = a(V) dt + p(V) dW. (4.2) 

If we explicitly modelled heat bath particles in Qq, then they would be distributed according to 
the Poisson distribution with density k^ in the interval (X(t) + R, oo). Their initial velocities would 
be given according to (2.2). Thus, using lemma 2.2 and (2.5), we can estimate the drift coefficent 
of the stochastic differential equation (4.2) to get 

2(v-V) H /X(t) + R- 
li + l V 



a(V) = 



At 



X(t)+R J 



At 



+ V - v) k^f^dvdx, 



where H(-) is the Heaviside step function. Using (2.2), we obtain 



a{V)-- 
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Va u V2 



exp 
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where k^ and are given by (2.1) and (2.2). In the limit jjl - 
we use the Taylor expansion in (4.3) to get 



(4.3) 

oo, we have V/^fjl + 1 -> 0. Thus, 



a(V): 



4V2 



/ny 



4j2D(n + 1) 



V z 



(4.4) 



The noise term in (4.2) can be computed by 

4(u - V) 2 



1 

At 



P 2 (V) 

Using (2.2), we obtain 
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Using (2.1), (2.2) and the Taylor expansion, we obtain 



Y 2 D+ * y^ y + 



3y 



2(ji + 1) 



-V 2 . 



(4.5) 
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Equations (4.4) and (4.5) are used in the multi-scale algorithm in table 3. The first two steps 
[Ml] and [M2] are the same as [Al] and [A2]. As heat bath particles are only simulated in the 
subdomain = (— L, 0), we remove all particles which left Q\) during the time interval (t, t + Af) 
in step [M3]. Step [M4] is the same as [A4], which introduces heat bath particles that have entered 
£2d through its left boundary x = —L during the time interval (t, t + At). The boundary at x = 0 is 
treated in step [M5] if the heavy particle does not intersect with this boundary. We assume that At 
is chosen so small that (2.7) is much smaller than 1. Then (2.7) can be interpreted as a probability 
of introducing one particle from the left (resp. right) during one time step. A particle introduced 
close to the right boundary of fi^ in step [M5] will have its position sampled according to the 
probability distribution 

for xe (-oo,0), (4.6) 



C\ erfc 



or M At^2 



where C\ is a normalization constant. The probability distribution (4.6) can be justified 
using the same argument as lemma 2.2 and equation (2.11). To sample from the probability 



Table 3. One iteration of the computer implementation of the multi-scale algorithm which is based on the MD model [A]. 



[M1] Compute 'free-flight positions' of heat bath particles and the heavy particle at time t + At using step [A1]. 

[M2] Compute post-collision velocities by (2.5) for every pair of particles which collided using step [A2]. 

[M3] Terminate trajectories of heat bath particles which left the subdomain £? D = (— L, 0). Update n accordingly. 

[M4] Implement the influx of heat bath particles through the boundary x = — L using step [A4]. 

[M5] If X(t) g (—R,R), then generate a random number r 2 uniformly distributed in (0,1). If r 2 < y 

(fi + l)Af/8, then increase n by 1, and introduce a new heat bath particle at position x n (t + At) with 
velocity v n (t + At) which are sampled according to probability distributions (4.6) and (4.7). 

[M6] If X(t) e (— R, R), then update the heavy particle velocity using (4.1). 

[M7] IfX(f) e [jR,L),then update the velocity of the heavy particle using (4.8). 

[M8] Continue with step [M1] using time t = t + At. 



distribution (4.6), we again use the acceptance-rejection algorithm in table 2 with parameters 
a\ and a 2 given by (2.10). In step [M5], we also sample the velocity v e R of the new particle using 
the truncated Gaussian distribution 



C 2 H(x-vAt)f ll (v), (4.7) 

where C 2 is a normalization constant. To sample random numbers according to the truncated 
normal distributions in steps [M4] and [M5], we again use the acceptance-rejection algorithm 
which is derived as proposition 2.3 in [18]. If the heavy particle does intersect with the boundary 
I, then step [M6] is executed. It uses (4.1) to incorporate collisions of heat bath particles from 
the right. If the particle does not intersect with Q^, then we simulate it in step [M7] using the 
discretized version of (1.7) given by 



V(t + At) = -yV(t)At + yV2DA££. (4.8) 

In figure 3, we present illustrative results computed by the algorithm [M1]-[M8]. We consider 
one heavy particle which starts at position X(0) = 0 with velocity V(0) = 0, as we did in figure 1. 
The distribution of its position at time t = l, computed using 10 5 realizations of the algorithm 
[M1]-[M8], is plotted in figure 3a. It is compared with the distribution obtained by the limiting 
BD model (2.6), which is, for t ^> y _1 , given by [19] 



y/47TD(t - t* 



exp 



x 2 



4D(f - t*) 



where t* = —. (4.9) 
2y 



In figure 3b, we plot the time evolution of the mean squared displacement. This figure can 
be directly compared with figure lb, because we use the same parameter values. The results 
computed by the multi-scale algorithm [M1]-[M8] compare well with the results given by the BD 
model (2.6). We have already shown in figure lb that the limiting BD model (2.6) also compares 
well with the MD simulations. In particular, the algorithm [M1]-[M8] is able to compute results 
with the MD-level precision by using coarser BD models in a part of the computational domain. 

The stochastic differential equation (4.2) was derived for collisions from the right. Using the 
same argument, we can also derive a stochastic differential equation which is approximating the 
effect of collisions from the left. We obtain 



dV = -a(-V) dt + p(-V) dW. (4.10) 

Adding (4.2) and (4.10) and using the independence of noise terms in (4.2) and (4.10), we can 
approximate collisions from both sides by the following stochastic differential equation for the 
velocity of the heavy particle: 



dV = (a(V) - a(-V)) dt + j P 2 (V) + P 2 (-V) dW. 



(4.11) 
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Figure 3. (a) Probability distribution of the heavy particle at time f = 1 computed by the multi-scale algorithm [M1]-[M8] 
(grey histogram) is compared with the distribution (4.9) given by the BD model (2.6) (black solid line), [b) The time evolution of 
the mean squared displacement computed by 10 5 realizations of the algorithm [M1HM8] (solid line) is compared with equation 
(2.13) (dashed line) and Vwi (dot-dashed line). We use /x = 10 3 , y =10,0 = 1, Af = 1(T 7 , L = 10, R = 1, X(Q) = 0 and 
|/(0) = 0. (Online version in colour.) 



Substituting (4.4) and (4.5), we derive (1.7). In particular, we have verified the limiting result in 
lemma 2.1. 



5. From three-dimensional moleculardynamics models [B] and [C] to Brownian 
dynamics 

We use a simple multi-scale geometry where domain Q = R 3 is divided into two half spaces. 
Heavy molecules are simulated in both half spaces. In = (— oo, 0) x R 2 , we use the MD model 
[B] or [C]. It is coupled with the BD model given by (1.6) and (1.7) in Qq = (0, oo) x R 2 . This set-up 
is a three-dimensional version of multi-scale problems which are schematically drawn in figure 2. 
Boundary conditions for heat bath particles at the interface I = {0} x R 2 can be specified using 
lemma 3.3. 

As in §4, we need to analyse the behaviour of a heavy molecule when it intersects with the 
interface I. Such a molecule is subject to the collisions with heat bath particles on the part of its 
surface which lies in fij^. This has to be compensated by using a suitable random force from Qq, so 
that the overall model is equivalent to (1.6) and (1.7) in the BD limit. To simplify the presentation 
of the algorithm, we use the same time step in and Qq. In §6, we present coupling of three- 
dimensional MD models with the BD model (1.1), which will make use of different time steps in 
different parts of the computational domain. 

The heavy particle is the ball with centre X = [Xi,X2,X3] with velocity V = [V\, Vi> ^3] aR d 
radius R. It intersects the interface I if X\(t) e (—R,R). Let us consider that heat bath particles are 
simulated in using the MD model [B] or the MD model [C]. Let us choose At so small that the 
probability of two collisions happening in the time interval (t, t + At) is negligible. As we do not 
explicitly simulate the heat bath particles in Qq, we will consider an additional correction of the 
velocity of the heavy particle in the form 

V(f + At) = V(f + At) + a(X(f), V(£)) At + P(X(t), V(t))VAt §, (5.1) 

where V(f + Af) is the post-collision velocity of the heavy particle at time t + At, which only takes 
into account collisions with the heat bath particles from £?d- It is either equal to V(£) or computed 
by (3.1) if a collision with a heat bath particle occurred in Q^. Note that we dropped the subscript 
\x in (3.1) and (3.2) to simplify our notation. Equation (5.1) is a generalization of (4.1) to three- 
dimensional simulations where a(X(f), V(t))At is the drift vector, P(X(t),V(t))y/~At % is the noise 



term and £ = is the vector of three normally distributed random numbers with zero 

mean and unit variances. Passing At -> 0, we observe that the contributions of the collisions from 
are given by the Ito stochastic differential equation 



dV = a(X(t), V(0) dt + P(X(t), V(t)) dW. 



(5.2) 



To estimate drift a and diffusion coefficient we separately consider MD models [B] and [C] 
in %ba and §5b, respectively. 



(a) Molecular dynamics model [B] 

The following lemma will be useful to estimate the drift coefficient a. 

Lemma 5.1. Let y > 0, D > 0,R> 0 and At > 0. Let us consider the MD model [B] where the positions 
and velocities of heat bath particles are distributed according to (3.3) and (3.4). Let us consider one heavy 
molecule in such a heat bath, i.e. N = 1, with the position of its centre to be at X(t) = [Xi(t),X2(t),X^(t)] 
and with velocity V(t) — \V\if), Viif), V^(t)]. Let y = (1/1,1/2,1/3) be a given point on the surface of the 
heavy molecule at time t, i.e. 



(yi - Xj(0) 2 + (y 2 - X 2 (t)f + (y 3 - X 3 (f)) 2 = R . 



(5.3) 



Then the average change of the j-th component of the velocity of the heavy molecule caused by collisions 
with heat bath particles in the time interval (t, t + At) at the surface area (y, y + dy) is iA)'(y) dy, where 



kpcrfaj-XjjfyAt | Ak^jyj-Xjjt^At 
(ji + 1)R 



V(f)-(y-X(f)) + 0(||V|| 2 ). (5.4) 



(fi + l)R 2 V2n 

Proof. Let us consider that a heat bath particle which was at point x at time t collided with the 
heavy molecule at time t + r e (t,t + At) at the surface point which had coordinate y at time t. 
Then the coordinate of the surface point at the collision time t + r was y + r V(f) and the pre- 
collision velocity of the heat bath molecule was v = V(f) + (y — x)/r. Using equation (3.1), we can 
write the change of the velocity of the heavy molecule during the collision as 



[v-V(f)] x = 



(y-x) (y-X(Q) \ (y-X(Q) 
+ t R ) R 



(5.5) 



The position x of the heat bath particle must be in the half space which lies above the plane tangent 
to the heavy molecule at the collision point y + rV(t). It can be parametrized by 



x = y + x V(f) + cir ^ + c 2 Trj2 + c 3 zrj 3 , 

where c\ > 0, C2 e R, C3 e R, and (y — X(t))/R, n\2, rj3 is the orthornormal basis in 
reads as follows: 



\ Then (5.5) 
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^2»?2 - ^3 ) dT dc 3 dc 2 dCi. 



(5.6) 



Substituting (3.4) for/ M and integrating over r, C2 and C3, we have 

kpty - Xj(t))AtV2 roo / vmn2 



^•(y) = - 



(/x + l)Ro> 



lci+V(Q ) exp 



c 2 

2^ 



dci. 



Integrating over ci, we deduce (5.4). 



Using lemma 5.1, we can compute the drift coefficient a(X(£), V(t)) in equation (5.2) as follows: 



aj(X(t),V(t))- 
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At 



s(x(0) 



(5.7) 



where S(X(t)) is the part of the surface of the heavy molecule which intersects the BD subdomain 
Qq, i.e. 

S(X(t)) = {y e £2c I y satisfies (5.3)}. 



Substituting (5.4) into (5.7), we have 



a/ (X(0,v(0) = - 



>-u.<y. 



(y ? -X(f))dy ,x 11 > y ' 

s(x(t))> 1 y (n + m 2 V2^ 



(H + 1)R 

Using (3.3) and (3.4) and evaluating the surface integrals, we obtain 



S(x(0) 



(y ; -X;(f)) 2 dy. 
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rYAA _ 3yA(/x + l)Dy / X?\ yVj I X\ 



and 



«;(X,V) = - 



2 + 3 



R R 3 



, for; = 2, 3, 



(5.8) 



(5.9) 



where we dropped the dependence on time t to shorten the resulting formulae. The noise matrix 
P(X(t), V(f)) will be estimated using /?(X(£),0), i.e. we will only use the first term in the Taylor 
expansion in V. Using similar arguments as in the proof of (5.6) and (5.7), we have 
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(5.10) 



for f = 1, 2, 3. Substituting (3.4) for/ M , (3.3) for A M and using /?(X, V) = /?(X, 0), we obtain 



X^ 

A,l(X,V) = yVD A /l + ^| 



and 



%(X,V) = yVD^/l + 
%(X,V) = 0, 



3Xi 



2R 2R 3 



, for; = 2,3 
for i =£j, 



(5.11) 



where the last equation can be verified using the same argument as equation (5.10). Note that, by 
substituting X\=R into (5.8), (5.9) and (5.11), we verify the limiting result in lemma 3.1. 



(b) Molecular dynamics model [C] 

Equations (5.6), (5.7) and (5.10), which are derived in %5a, are applicable to both MD models [B] 
and [C]. To estimate the drift coefficient a(X, V) for the MD model [C], we substitute (3.5) for A M 
and (3.6) for/ M in (5.6) and (5.7). We obtain 



and g?2(X,V) and a^(X,V) are again given by (5.9). Substituting (3.5) and (3.6) in (5.10) and 
integrating, we obtain that noise matrix /?(X, V) satisfies (5.11). We again note that the special 
choice Xi = R in (5.12) can be used to verify the limiting result in lemma 3.2. 
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Figure 4. (a) Probability distribution of the first coordinate, X h of the heavy particle at time t = 1 computed by the three- 
dimensional multi-scale algorithm (grey histogram) is compared with the distribution (4.9) given by the BD model (black solid 
line), (b) The time evolution of the mean squared displacement computed by 10 5 realizations of the three-dimensional multi- 
scale algorithm (solid line) is compared with the limiting BD model (1.6) and (1.7) (dashed line) and -JWt (dot-dashed line). 
We use ii = 10 3 , y = 10, D = 1, At = 10" 6 , L = 5, R = 1, X(0) = [0, 0, 0] and V(0) = [0, 0, 0]. (Online version in colour.) 



(c) Illustrative numerical results 

In %ba,b, we have observed that the only difference between MD models [B] and [C] is a different 
formula for the coefficient «i(X, V) in (5.2), given by (5.8) and (5.12), respectively. The remaining 
terms in (5.2) are the same, given by (5.9) and (5.11). In this section, we present an illustrative 
computation with the MD model [C], but the same results can also be obtained with the MD 
model [B] (results not shown). An illustrative computation with the MD model [B] is presented 
later in §6. 

We consider a three-dimensional generalization of the illustrative problem from figure 3. One 
heavy particle which starts at position X(0) = [0, 0, 0] with velocity V(0) = [0, 0, 0] is simulated 
using a three-dimensional generalization of the algorithm [M1]-[M8]. We use the MD model [C] 
in ^ D = (-oo,0) x R 2 and the BD model (1.6) and (1.7) in Q c = (0,oo) x R 2 . In step [M6], we 
replace (4.1) with its three-dimensional analogue (5.1), where drift a and diffusion coefficient 
P are given by (5.12), (5.9) and (5.11). The distribution of Xi positions of the heavy particle at 
time t = l, computed using 10 5 realizations of the multi-scale algorithm, is plotted in figure 4#. 
The limiting BD result is again given by (4.9). In figure 4b, we plot the time evolution of the 
mean squared displacement. The mean squared displacement corresponding to the limiting BD 
model (1.6) and (1.7) is given by (2.13) multiplied by V3 because we have three spatial dimensions 
(dashed line). As expected, the models compare well. The mean squared displacement obtained 
for the BD model (1.1) is plotted as the dot-dashed line. 



6. Application to protein binding to receptors 

In this section, we apply our results to a simplified model of protein binding to receptors on 
the cell membrane. We consider simple geometry, which is schematically shown in figure 5. Our 
computational domain is a part of the intracellular space next to the cell membrane given as the 
cuboid Q = [0,Li] x [0,L2] x [0,L2], where L\ > 0 and L2 > 0. The cell membrane is modelled by 
one side of the cuboid, namely 

3^M = {0} x [0,L 2 ] x [0,L 2 ], 

which is shaded grey in figure 5. Our goal is to model the binding of diffusing proteins to 
receptors on the cell membrane with an MD level of detail. Therefore, we define £2d as a part 
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Figure 5. Schematic of the computational domain used in the protein binding example. (Online version in colour.) 



of the intracellular space which is close to the cell membrane dfiyi, i.e. 

^ D = [0, h] x [0, L 2 ] x [0, L 2 ] and Q c = IK L\] x [0, L 2 ] x [0, L 2 ], 

where h > 0 and the interface J is at x\ — h. Diffusing proteins are modelled as spheres of radius R. 
We consider that a protein which hits the boundary dQyi will bind to a receptor with probability 
P, and otherwise it is reflected. This type of a reactive boundary condition is common for BD 
simulations [20]. In the case of MD, more detailed models of protein binding could be introduced 
in &v [21,22]. However, the main goal of this section is to show how an MD model in fi^ 
can be coupled with BD simulators which have been developed for simulations of intracellular 
processes. Therefore, we keep the MD model in as simple as possible. 

If we used BD model (1.6)-(1.7) in Qq, then the situation would be more or less the same 
as in §5. However, modern BD simulators of intracellular processes work with the high-friction 
limit (1.1) rather than (1.6)-(1.7). For example, the software package Smoldyn discretizes (1.1) 
with a fixed time step and uses (1.2) to update positions of diffusing proteins. In particular, it uses 
larger values of time step than we used in §5. Then the problem can be formulated as follows: we 
would like to use the MD model with time step At in Q\) and couple it with the BD model (1.2) 
with larger time step At, namely 

Xi(t + At) = Xi(t) + y/lDAt^, i = 1,2,3, (6.1) 

if the diffusing molecule is far away from Q^. We couple these models using the intermediate BD 
model (1.6)-(1.7). We introduce two additional interfaces 

h = {hi) x [0, L 2 ] x [0, L 2 ] and I 3 = {h 3 } x [0, L 2 ] x [0, L 2 ], 

where h < /z 2 < /13 < L\, as shown in figure 5. We denote 

Q C1 = [h,h 3 ] x [0,L 2 ] x [0,L 2 ] and Q C 2 = \hiM\ x [0,L 2 ] x [0,L 2 ], 

i.e. Qq\ and ^?C2 are two overlapping subdomains of Qq. We simulate the time evolution of 
the position X(f) of one protein molecule. If X(£) e Qq2, then the BD model (6.1) will be used in 
Qq2 until the molecule leaves Qq2- Then we switch to the shorter time step At and use the BD 
model (1.6)-(1.7) in Qq\. The protein molecule can leave Qq\ in two possible ways: 

(i) The protein molecule crosses the interface I3. 

Then we revert to the BD model (6.1) which is used in £?C2- 

(ii) The protein molecule crosses the interface I. 

Then we use the method from §5 for coupling the MD model in with the BD 
model (1.6)-(1.7) in Q C1 . 

As the subdomains Qqi and Qqi overlap, we can use the limiting result (4.9), which implies 
that, for times t~>y~ x , the BD model (1.6)-(1.7) is given by the BD model (6.1) shifted by time 
t* =3/(2y). In particular, we will also add or subtract t* from the time variable whenever we 
switch between BD models. 




Figure 6. (a) Distribution of positions along the X r axis at time t = 1 computed by the multi-scale model described in §6 (grey 
histogram). The black solid line is the solution of the limiting PDE model (6.2) and (6.3). Vertical lines denote interfaces /, l 2 and 
/ 3 . [b) Probability that the protein is bound to a receptor as a function of time t computed by the multi-scale model (solid line) 
and the PDE model (6.2) and (6.3) (dashed line). Parameters used: D = 10, y = 10 2 , \i = 10 3 , K = 1, R = 1, k = 10 2 , h = 4, 
h 2 = 8, /? 3 = 12 and /? 4 = 3. (Online version in colour.) 



In figure 6, we present illustrative results computed by averaging over 10 5 realizations. Initial 
positions of the protein molecule are uniformly distributed along the Xi-axis. The histogram of 
positions (along the Xi-axis) at time t = 1 is plotted in figure 6a. Interfaces I, I2 and J 3 are also 
shown in this plot. As we used a very simple model of the protein binding, we can compare it 
with the mean-field limit given by the solution of the PDE 

with boundary conditions [20] 



(x 1 ,t) = D T ±(x 1 ,t), x 1 e[0,L 1 l t>0, (6.2) 



D — (0,t) = KQ(0,t) and D—(L 1 ,t) = 0, whereP=^2^. (6.3) 
dxi dxi ^/Dy 

The solution of (6.2) and (6.3) with uniform initial condition q(x\, 0) = const, is given by the 
black solid line in figure 6a. As we only visualize the distribution along the Xi-axis in figure 6a, 
we can further decrease the computational cost by truncating the simulation domain in the x-i 
and X3 directions to the region close to the protein molecule. That is, we only simulate small 
particles in the subdomain [Q,h] x [X 2 (f ) - h±,X 2 (t f ) + h A ] x [X 3 (f) - /i4,X 3 (f) + /14], where I' is 
the time when the protein molecule enters £?d u ^ci • This subdomain (moving window) is 
shifted accordingly whenever X2H) or X 3 (f) approach its boundary. 

The probability that the protein is adsorbed to the surface is given as a function of time in 
figure 6b. It again compares well with the results obtained by the limiting PDE system (6.2) 
and (6.3). 



7. Discussion 

In this paper, a multi-scale approach that uses MD simulations in a part of the computational 
domain and BD simulations in the rest of the domain has been presented and analysed. The 
ultimate goal of this research is to use MD to help parametrize BD models of intracellular 
processes. One application area is modelling proteins in an aquatic environment, which is useful 
for understanding protein binding to receptors (surfaces) as shown in §6. 

The main idea of the presented coupling of MD and BD models is based on using 
equations (4.1) and (5.1) and estimating drift and diffusion coefficients for velocities of molecules 
which cross the interface I. This coupling uses the same time step for the BD model (1.6)-(1.7) as 



for the MD model. In §6, it was shown that this is not a limiting step of this approach, because the 
BD model (1.6)-(1.7) is only needed in a small part of the domain next to £2^. Then the coarser 
BD model (6.1) with larger time step can be used in the rest of the simulation domain, using a 
suitable overlap region. Another overlap region could be used to couple BD simulations with 
mean-field PDE-based models [11]. Then multi-scale models which couple BD (of point particles) 
with coarser reaction-diffusion approaches would be capable of further increasing time scales 
and space scales of simulations [10,11]. 

MD models considered in this paper are relatively simple and analytically tractable, describing 
water molecules as point particles. An important generalization is to consider more complicated 
MD models of water molecules [23]. For example, Rahman & Stillinger [24] model water 
molecules as rigid asymmetric rotors. That is, each water molecule is described by six coordinates: 
the position of its centre of mass and three angles describing molecule orientation. The energy 
of water solution is given as the sum of kinetic energies (for translation and rotation) and the 
intermolecular potential which is assumed to be pairwise additive and can be given in several 
different ways, i.e. the heat bath is given by its Hamiltonian [23-25]. I am currently investigating 
MD models based on Hamiltonian dynamics, with the aim of designing and analysing multi-scale 
algorithms similar to the algorithm [M1]-[M8] from this paper. The ultimate goal of this research 
is to design BD models of intracellular processes which make use of modern MD simulations 
[26,27] in small regions of intracellular space where a BD description is not available or applicable 
(and an MD level of detail is required). Results will be reported in a future publication. 
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